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Abstract: We describe an algorithm for long-term planetary orbit in- 
tegrations, including the dominant post-Newtonian effects, that employs 
individual timesteps for each planet. The algorithm is symplectic and ex- 
hibits short-term errors that are 0(eSl 2 r 2 ) where r is the timestep, £1 is 
a typical orbital frequency, and e <C 1 is a typical planetary mass in solar 
units. By a special starting procedure long-term errors over an integration 
interval T can be reduced to 0(e 2 Sl 3 r 2 T). A sample 0.8 Myr integration of 
the nine planets illustrates that Pluto can have a timestep more than 100 
times Mercury's, without dominating the positional error. Our algorithm 
is applicable to other A-body systems. 
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1. INTRODUCTION 

The longest integrations of planetary orbits are still well short of the age 
of the solar system. So far the full planetary system has been followed for 
100 Myr (Sussman & Wisdom 1992), and the five outer planets for 1 Gyr 
(Wisdom & Holman 1991; hereafter WH). A semi-analytic secular pertur- 
bation theory has been used to follow the planetary system for 200 Myr 
(Laskar 1989, 1990) and shows very good agreement with direct integra- 
tions (Laskar et al. 1992, Sussman & Wisdom 1992). Progressively longer 
integrations have generally revealed interesting new phenomena, notably 
weak chaos in the orbits of the inner (Laskar 1989, 1990) and outer (Suss- 
mann & Wisdom 1992) planets. This situation motivates the development 
of faster and more accurate integration methods. 

Traditionally, long solar system integrations have used high-order mul- 
tistep integration methods in Cartesian coordinates (see, for example, 
Quinn et al. 1991, hereafter QTD). However, substantial improvements in 
speed are possible using integration methods that are specifically designed 
for motion that is (a) Hamiltonian, and (b) nearly Keplerian; we call these 
mixed variable symplectic integrators and they are the subject of this paper. 

The mixed variable symplectic (or MVS) integrators were introduced 
by WH, and also (independently and by different arguments) by Kinoshita 
et al. (1991; hereafter KYN). These derive their advantage by switching con- 
tinually between Cartesian variables (wherein the perturbations are easy to 
evaluate) and Kepler elements (which make the solar part simple) — hence 
'mixed variable'. The symplectic property (i.e., having certain Hamiltonian 
conservation laws built-in) helps control long-term errors. Saha & Tremaine 
(1992; hereafter Paper I) describe a startup technique ('warmup') that sub- 
stantially improves the long-term accuracy of MVS integrators; a related 
technique ('symplectic correctors') is given by Wisdom et al. (1994; here- 
after WHT). In this paper we generally follow KYN's methods of analysis 
but WH's algorithms. 

A limitation of symplectic integrators so far is that they generally 
allow neither adaptive nor individual timesteps. Adaptive timesteps are 
indispensable for situations with close encounters or very high eccentricities, 
but for planetary and satellite orbits we do not miss them much. However, 
the requirement that all planets be followed with a common timestep is 
certainly undesirable, since planetary orbital periods range over three orders 
of magnitude. If the common timestep is dictated by the desired accuracy 
for Mercury, then for Pluto we may be paying for several unnecessary orders 
of magnitude of accuracy (see Fig. 2). 
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The main contribution of this paper is to introduce an MVS integrator 
with individual timesteps. We describe our new algorithm in Sec. 3 after 
covering some operator formalism in Sec. 2. Section 4 has details of the 
equations of motion. In Sec. 5 we show how the leading order general 
relativistic corrections can be neatly incorporated. Section 6 has some 
numerical tests. Finally, in Sec. 7 we discuss some variations on our method. 

The savings in computer time from adopting individual timesteps is 
modest in the case of the planetary system (a little over a factor of two), 
but still significant considering that solar system integrations often run for 
weeks or months of machine time. For other problems, however, (e.g. if the 
lunar orbit is integrated as well) the savings can be much larger. 



2. LEAPFROG 

Hamilton's equations can be written in terms of a Poisson bracket operator 
as 

which has the formal solution 

e*< (2) 

An integration algorithm can be thought of as an approximate expression 
for the operator in (2). 

A common variety of Hamiltonian is the sum of two (or more) parts, 
each of which is soluble in isolation. That is, 

H = H A + H B , (3) 

where ,Ha ^ and ,Hb ^ are known. These two operators can be used as 
building blocks to construct approximations to (2). The obvious example 
is 

2 

H A =^, H B = V(r), (4) 

as a kinetic energy and coordinate-dependent potential. 

The solar system Hamiltonian can be expressed in the form (4), but it 
is more useful to take 



Ha — Hjiep, Hb — Hi D t 



(5) 
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where i?Kep is the kinetic energy plus solar potential and _ff; nt is the inter- 
action energy between the planets (WH, KYN). Then e*^ ,f?Ktp ' generates 
motion along unperturbed Kepler orbits, while e*^ > H ^} generates a change 
of momenta with the coordinates fixed. 

A key ingredient in integration algorithms for such systems is the op- 
erator identity (see, for example, Yoshida 1993): 

e \r{ ,H A ] e r{ ,H B } e \r{ ,H A ] = g r{ ,ff A+ ff B+ ff err } ^ 

where H eTT is a formal power series in r starting at 0(t 2 ) and consisting of 
nested Poisson brackets of Ha and Hb'- 

Herr = ^{{#A, #B } , #B + \H A } + 0(r 4 ). (7) 

In general the series for H eTT does not converge 3 , and is interpreted 
as an asymptotic series. The left side of (6) is recognized as one step of a 
generalized form of leapfrog, with r being the timestep. The well-known 
result that leapfrog is second order follows from the fact that H eTT is 0(r 2 ). 
The expression (6) also reveals other useful properties: 

(i) The integration errors are Hamiltonian; that is, the integration algo- 
rithm follows exactly the dynamics of a nearby 'surrogate' Hamiltonian 

H + H eTT . (8) 

This result must be interpreted cautiously since H eTT is only a formal 
series, but in practice, for small r, analyses based on the leading term 
in if err provide considerable insight. 

(ii) For orbital motion, (i) suggests that the energy error is bounded and 
therefore the position errors grow as the integration time T, rather 
than as T 2 (as in most integration methods). Moreover, the rate of 
error growth can be estimated from H eTT , as illustrated by KYN. 

Now suppose that Hb is 0(e) smaller than Ha, as in planetary motion if 
we identify Ha and Hb as in Eq. (5) (in this case e is of order the planetary 
mass). Two further properties then follow. 



3 This is illustrated by the following example (J. Wisdom, private communication): 
consider the pendulum Hamiltonian, with Ha = ^p 2 and = — cosg. The map 
(6) applied to z = (q,p) is then simply the standard map (to within a shear). The 
standard map exhibits chaos, which is incompatible with motion in an autonomous one 
degree-of- freedom Hamiltonian Ha + H& + H eTT . 
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(iii) Since H eTT consists entirely of nested Poisson brackets of Ha and Hb , 
H eTT will be of 0(er 2 ). 

(iv) The long term errors can be reduced from 0(ct 2 T) to 0(e 2 T 2 T) by 
special starting procedures (Paper I — see also WHT). The reason is that 
if err has no secular terms at 0(e), 4 so that to this order the relation 
between actions and frequencies is the same in the actual and surrogate 
systems. Thus any 0(er n T) error terms come from a difference in 
the values for the actions in the actual and surrogate systems — this 
difference of course giving rise to a constant frequency error and hence a 
linearly growing position error. If the difference between the surrogate 
and exact action values can be removed (to leading order in e) by a 
suitable small alteration in the initial conditions, all 0(er n T) error can 
be suppressed. 

By concatenating leapfrog steps one can produce higher order integra- 
tors; for example, three consecutive leapfrog steps with timesteps in the 
ratio 1 : — 2~ : 1 amount to a single step of a fourth order integrator. 
Yoshida (1990) shows that arbitrarily high orders are possible. With suit- 
able modifications, properties (i)-(iv) carry over to higher order integrators. 

In this paper, we develop an integration algorithm with individual 
timesteps for each planet. The idea is to apply Eq. (6) recursively, replacing 
the operator e*^ ,Hb ^ by a more complicated operator that itself involves a 
leapfrog step. The arguments leading to properties (i)-(iv) are not affected 
by this change. 



3. LEAPFROG WITH INDIVIDUAL TIMESTEPS 

Suppose that the Hamiltonian for a system of planets is split into Hj^ ep and 
Hi n t as in Eq. (5). The details of what is inside Hj^ ep and Hint we leave 
for Sec. 4. Now imagine two clocks K and I, associated with Hj^ ep and 
Hi n t respectively. These resemble the clocks used in chess tournaments in 
that only one of them is running at a given time; when the 'Kepler clock' 
K is running, each of the planets moves along its osculating Kepler orbit, 
and the interplanetary interactions are turned off; when the 'interaction 
clock' I is running all the coordinates stay fixed while the momenta change 



4 The exposition in Sec. 3 of Paper I is flawed in that it shows only that the error 
Hamiltonian has no secular terms at 0(er 2 ) (P.-V. Koseleff, private communication); 
however, the argument in that Section readily extends to show that there are no secular 
terms of order 0(er n ) for any n. 
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according to H ln t- Thus an interval r of A' or I corresponds to the operators 
g r{ ,-ffRep} or g T{ ,H int } reS p ec tively. A single leapfrog step of length r can 
be written in pseudocode as the following procedure: 

(Advance K by ^r) 

( Advance I by r ) (9) 
(Advance K by ^r) 

We can write a sequence of leapfrog steps, each of size r, in terms of the K 
and I clocks as follows. 

(Advance K by ^r) 
loop 

( Advance I by r ) 

( If output is desired then exit loop) (10) 
( Advance K by r ) 
end loop 

(Advance K by -|t) 

Note that K and I show the same time at the start and end, but in between 
they are not synchronized. 

Now we propose a generalization to individual timesteps. First we 
assign each planet its own timestep (which normally does not vary during 
the integration). Assume that the planets are indexed from innermost 
outwards as 1. .N , planet i has timestep r;, and we restrict t;+i to be an 
integer multiple of r;. We assume that Hj^ ep and H ln t can be written in 
the form 

N N 
#Kep = ^ -ffKep, i, H ; n t = ^ H ; n t , i ■ (11) 

8 = 1 8 = 1 

Here -ffxep, i is the Kepler Hamiltonian for the two-body system consisting 
of the sun and planet i, and H ln t t i is the potential energy arising from the 
interaction of planet i with planets i + 1 to N . These Hamiltonians have 
the properties 

{-ffKep.i, #Kep,j} = 0, {-ffint.8, -ffint.j} = 0, (12a) 



and 



{H Ke p,i,H mt j} = 0, for j > i. 



(12b) 
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(In Sec. 7 we discuss algorithms for the case where Eq. 12b does not apply.) 

We can now assign clocks 7Q and 7 8 - to each planet: advancing a Kepler 
clock Ki by r corresponds to the operator e T>l - ,J?Kep ''', that is, to advancing 
planet i along its osculating Kepler orbit; advancing an interaction clock 
Ii corresponds to the operator e T>l - ' Hint - that is, to changing momenta 
according to the interactions between planet i and all planets j > i (note 
that the clock 7jv does nothing). Equations (12a, b) imply that the order of 
advancing any two Kepler clocks Ki, Kj can be reversed without affecting 
the result (i.e. the corresponding operators commute); the same is true for 
any two interaction clocks Ii, Ij, and also for Ki, Ij when i < j. 

The general idea of our algorithm is to interleave advances of the vari- 
ous clocks such that (i) at the end of the integration each of the clocks has 
been advanced by the same amount; (ii) the integration is reversible, that 
is, if the algorithm is applied to the time-reversed final state we recover 
the time-reversed initial state; (iii) the clocks remain as close to synchro- 
nization as possible. A suitable algorithm is expressed by the following 
recursive procedure: 

procedure TICK (i) 

(Advance Ki by ) 
( Advance Ii by r; ) 
if i > 1 

loop r 8 /r 8 _i times 

call TICK (i - 1) 
end loop 
end if 

(Advance K{ by ) 
end TICK 

The order of the step that advances Ii and the loop that calls TICK (i — 1) 
can be reversed without affecting the result; bearing this in mind it is easy 
to see that the algorithm is time-reversible. To advance the integration of 
the TV-body system by rjv one simply calls TICK (N). 

However (13) is not the most efficient form of the algorithm, because 
it often splits an advance of a Ki clock by r; into two successive advances 
by \ti. The form we actually implement is an equivalent non-recursive 



(13) 
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version: 

(Advance all the K{ by ) 
loop 

( For any i where Ki has changed more recently than Ii , 
advance Ii by r; ) 

( If the Ii are all equal and output is desired then exit 
loop) ( 14 ) 

loop for i = 1. .N 

( If i = 1, or Ki + \ti < Ki-i, advance Ki by r; ) 
end loop 
end loop 

(Advance all the Ki by ) 

Note that although the clock 7jv does nothing, it must be monitored along 
with the other interaction clocks. It may be helpful to follow the steps in 
the algorithm by hand in a few simple cases to see how it works. 

As an example, consider the case of two planets, with = 2ti = 2r. 
Since I2 does nothing we may write I\ = I. A single step (of duration 2r) 
is: 

(Advance K\ by ±r and A'2 by r) 
( Advance I by r ) 
( Advance K\ by r ) 
( Advance I by r ) 
(Advance K\ by and A'2 by r) 

It is straightforward to derive the error Hamiltonian for this case: 

H eTT = ^ [(17, 7) + ±(17, 1) + 4(27, 7) + 2(27, 2) + 4(27, 1)] + 0(r 4 ), 

(16) 

where the symbols (ij, k) represent nested Poisson brackets; for example, 
(27, 1) = {{77Kep, 2, -ffint}, -ffKep, l}- The first two terms in Eq. (16) repre- 
sent the error arising from the leapfrog step of length r in planet 1 (cf. Eq. 
7); the second two arise from the leapfrog step of length 2r in planet 2; and 
the final term is associated with the presence of both planets. 

Notice that when the interaction between two planets i and j is com- 
puted, the clocks Ki and Kj are generally at different epochs. Second-order 



(15) 
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accuracy is still achieved because on average K{ is behind Kj as much as 
it is ahead of Kj . Accuracy could be improved by synchronizing the K{ 
clocks when the interactions are to be computed. One way to do this would 
be to run all the Kj^i clocks until they synchronize with K{, advance J; by 
T{, then run all the Kj^i clocks back to their previous settings. This strat- 
egy is not useful for solar system integrations, because moving the planets 
back and forth along Kepler orbits takes too much computer time to be 
worth the gain in accuracy (it would be useful in a system where advancing 
K{ took much less computer time than advancing /;). A better plan is to 
replace the Kepler orbit, for the purpose of the synchronization only, by 
an approximation (we might call this 'symplectic interpolation'). Even a 
crude approximation, provided it is symplectic, improves the accuracy sub- 
stantially. We have found that a simple rotation in the invariable plane by 
a preset amount works well. That is, we replace advances of any J; clock 
with the sequence of steps: 

(Evaluate <f>j = iij(Ki — Kj) for j > i, hj being the mean 
motion of planet j at the start of the integration ) 

( Rotate each planet j > i in the invariable plane by <f>j ) (17) 

( Advance J; by r; ) 

( Rotate each planet j > i in the invariable plane by —<f>j ) 



As described so far the algorithm leaves errors of 0(er 2 T). To reduce 
long-term errors further to 0(e 2 T 2 T) we use a special starting procedure. 
The general idea is to change adiabatically (i.e., over a time much longer 
than the orbital time but much shorter than the total integration time) from 
a much more accurate integration procedure to the one we will actually 
use. This procedure ensures that the error Hamiltonian H eTT grows slowly 
so that the actions are unchanged (cf. point (iv) in Sec. 2), which removes 
the leading source of long-term error. In Paper I we recommended starting 
with a very small timestep and then gradually increasing the timestep to 
its final value. (With individual timesteps, this would require starting with 
all the T{ scaled down and then gradually scaling them up again, always 
maintaining a fixed ratio between the r 8 .) This procedure did not produce 
the hoped-for gain in accuracy with our new algorithm — apparently be- 
cause changing the timestep causes the system to sweep through artificial 
resonances (see Wisdom & Holman 1992), which are prominent because of 
the large timesteps used by some of the planets. A slightly more subtle 
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startup procedure works better: 

( Integrate backward (or forward) for a large number of orbits 
with the T{ all scaled to small values (by at least ~ s/^), while 
gradually reducing the interaction strengths to zero ) 

( Integrate forward (or backward) to the starting point with 
the regular timesteps r;, while gradually reviving the inter- 
action strengths to the correct values ) 

During the first part H eTT is negligible (because the timesteps are small); 
during the second part H eTT is initially negligible (because the interplane- 
tary interactions are small and the integrator follows Kepler orbits exactly) 
and then grows slowly, as required for adiabatic invariance of the actions. 
Hereafter, we refer to this procedure as a 'warm start' or 'warmup'. Its 
computational overhead is small compared to the main integration. 

The warmup technique is based on eliminating any 0(e) contribution 
from iferr to the actions. The symplectic corrector approach of WHT is 
to annihilate the 0(e) part of H eTT altogether with a canonical transfor- 
mation. WHT apply the transformation to the initial conditions, integrate 
(using equal timesteps) as usual, and then apply the inverse transformation 
whenever output is desired. The difference between symplectic correctors 
and warmup is analogous to the difference between perturbation theory and 
averaging. As one might expect, both methods are equally good at control- 
ling long-term errors; symplectic correctors have the advantages that they 
also remove short-term oscillatory errors at 0(e), and that they can be 
extended to higher order in e. They are, however, more complicated to im- 
plement, especially if extended to higher orders or to individual timesteps. 



(18) 



4. EQUATIONS OF MOTION 

The material in this section mostly follows WH's discussion; however, we 
include a little more detail on the calculation of time-evolution under if; n t, 
and skip over some other points not essential in our context. 

The MVS integrators expect a Hamiltonian expressed as a sum of Ke- 
pler terms of the type p 2 /2m — fi/r and interaction terms of the type V(r); 
moreover, for the integrator to work efficiently the interaction terms asso- 
ciated with each body must be much smaller than its Kepler term. These 
requirements necessitate a special set of variables (well known in classical 
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perturbation theory and for the same reasons), the Jacobi variables. Helio- 
centric variables will not do because then the Hamiltonian does not have the 
right form, and barycentric variables will not do because the sun's motion 
is not dominated by a Kepler part. To convert to Jacobi coordinates one 
orders the planets (inner-to-outer being usual and probably best, but not 
essential) and reckons the coordinates of each planet from the barycenter 
of the sun and all the previous ones. 

We use Gaussian units: the solar mass is unity and k 2 is the gravita- 
tional constant. The planets have masses m; and heliocentric coordinates 
r; (in this section the dummy indices i and j always range from 1 to N). 
We first define cumulative masses o - ;, and renormalized masses and gravi- 
tational constants m; and nf. 



<Ti = (Ti-i +rrii, (T = 1, 



m; = m; 



(19) 



The Jacobi coordinates f; are then defined by 



r,; = r, : 



which has the 



\ - rnjYj 



j<i 



<7,' 



(20b) 



If we add to the set fi. .fjv the position (fo, say) of the barycenter of the 
whole system (in some inertial frame), then the momenta (po. -Pn, say) 
conjugate to fo. .fjv have the simple interpretation that po is the total 
momentum and ^ 

Pi = fhiVi, where v; = . (21) 

The canonical set fi. .fjv,Pi- Pn is collectively known as the Jacobi vari- 
ables; fo is ignorable and we disregard fo,Po from now on. 

The advantage of transforming to Jacobi variables is that in the 
barycentric frame 

~2 

Kinetic energy = — ^— . (22) 
^-^ 2 m,- 
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For a derivation see Plummer (1960). WH interpret the transformation in 
an interesting way as a matrix diagonalization. In view of (21), we can 
write the full Hamiltonian as 



where (cf. Eq. 11) 



H — i?Kep + -ffint + #misc (23a) 



~2 ^ 

-ffKep = / Hj£ep,i , where i?Kep,i = 77-! Vi^—, 

■^-^ Znii ri 

H m t = ^ -ffint, i , where H mtt ; = i?dir, i + -ffindir ^'1 , 

i 



(23b) 



r,- — r,- 

i>i 1 J 



-ffmdir = k 2 V" TTli I ^- - — 

' J \ J' ' J' ' 



The designations i?dir and -ff; n dir refer to the direct and indirect parts in 
the usage of WH; this is similar but not identical to the traditional usage. 
-Hmisc denotes any other physical effects. The most interesting of these is 
the general relativistic correction, which we will approximate with a post- 
Newtonian Hamiltonian ifpN to be discussed in the next section. Besides 
this we include a term H\ un that represents the attraction of the sun on the 
quadrupole moment of the Earth-Moon system (see QTD for details). Thus 
we have 

-ffmisc = -ffpN + -fflun- (23c) 

Other effects such as solar oblateness, asteroids, and the galactic tidal field 
are thought to be < 10 -10 of the Kepler part (see QTD), and we neglect 
them. 

Now we consider time evolution under i?Kep and _ff; nt . 

i?Kep of course generates evolution along a Kepler orbit. Computing 
this involves (implicitly or explicitly) transforming from Cartesian positions 
and velocities to orbital elements, incrementing the mean anomaly by the 
appropriate amount, and then changing back to Cartesian variables. As 
advocated by WH, this process is efficiently encapsulated in Gauss's / and 
g functions. 
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Now consider ifdir- It is convenient first to compute the impulse in 
heliocentric variables and then to transform to Jacobi variables. The accel- 
erations due to ifdir i are 



dt 
dt 



dir. i dt Idir, i — l 



-E- 



dvj_ 

dt 



(24) 



dir. i 



Evolution under i?; n dir is more complicated. Because i?; n dir mixes r; 
and ?i, it cannot be usefully be split into a sum of contributions associated 
with each planet. The acceleration due to i?; n dir is 



dvj 
dt 



indir 



l 



(25) 



We lump i?indir in with i?; n t, l for the algorithm, and thus compute Eq. (25) 
whenever the I\ clock is advanced. This procedure does not produce a 
bottleneck, because all the dv i / dt[ lnt iir can be computed in O(N) opera- 
tions, while the d\i/dt\ c \\ Tt ; computations require 0(N 2 ) operations. We 
should also mention the standard practice of rewriting differences such as 
ii/rf — Ti/rf in (25) to be less sensitive to roundoff error — see, for example, 
the discussion of Encke's method in Danby (1988). 



5. POST-NEWTONIAN CORRECTIONS 

General relativistic effects in planetary motion have fractional amplitude of 
order k 2 /c 2 r ~ 10 -8 at r = lAU. Neglecting corrections that are smaller 
still by 0{k 2 / c 2 r) or by O(m), the relativistic effects can be expressed 
through the post-Newtonian Hamiltonian (see Landau k Lifshitz 1975) 5 : 

H ™ --iL.y-^f ~ g^f " > (26) 



5 This form for //pn assumes that the metric has the form ds 2 = [1 — 2k 2 / (c 2 r)]dt 2 — 
[1 + 2k 2 1 (c 2 r)]rfx 2 + 0(c -4 ), consistent with the isotropic or harmonic form of the 
Schwarzschild metric. This is the metric recommended by the IAU, but older papers 
(e.g. Brouwer & Clemence 1961) often use the "standard" form of the Schwarzschild 
metric, which yields different equations of motion. 
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note that the distinction between barycentric, heliocentric, and Jacobi co- 
ordinates, and between fi{ and k 2 or m; and m;, is negligible to the accuracy 
we are considering. 

ifpN mixes coordinate and momentum dependencies and hence is not 
easily decomposed in the form (3). The usual practice has been to re- 
place it by an ad hoc potential Upn designed to mimic the most important 
effects (Nobili k Roxburgh 1986, Paper I, Sussman k Wisdom 1992). Al- 
ternatively one could replace ifpN by its average along Kepler orbits ifpN- 
Neither of these approaches is entirely satisfactory. The reason is that 
the difference between ifpN and ?7pn or ifpN usually contains a secular 
part; this does not noticeably affect the perihelion precession, but gives the 
orbital frequency an error of order the post-Newtonian effect itself (this 
could perhaps be alleviated by some special starting procedure analogous 
to warmup, though it is not obvious how). 

However there is a simple extension of the MVS method that easily 
accommodates ifpN, and to which we devote the rest of this section. We 
rearrange Eq. (26) as 

#pn = J2 ( a ^Kep, i + Pi/rj + TiPi ) (27a) 

i 

where 

ai = 3/(2m 8 -c 2 ), p i = -n 2 /c 2 , 7i = -l/(2mfc 2 ). (27b) 

This is similar to the form (3) — except that there are now three terms, 
each of which is integrable in isolation — and we may thus compute time 
evolution as follows. 

Each tti-ffxep i we combine with the corresponding i?Kep, i in Eq. (23b). 
Then advancing the K{ clock by r; is redefined as the operator 

exp (t;{ , if K ep, i + ai^Kep, J) • (28) 

Under the operator (28) H^ ep i is conserved and equals —^/lifhi/ai (a; be- 
ing the osculating semi-major axis in Jacobi variables), so (28) is equivalent 
to 

exp (r/{ ,H Kep>i }), ri=(l- ^) r„ (29) 

Thus, to incorporate the aiH^ i terms, we merely rescale the time argu- 
ment passed to the / and g functions. The Pi/f 2 are trivial to deal with — we 
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The jipf terms we incorporate through leapfrog 
operators that evolve under a pf term before and after advancing a A' s - 
clock. 

Because ifpN changes the momentum dependence of the total Hamil- 
tonian, we no longer have p 8 = fhidvi/dt (Eq. 21), so the expressions for 
accelerations in Sec. 4 have to be modified. A simple way to incorporate 
the necessary modifications is as follows. We redefine v 8 - to be a pseudo 
velocity: 

(ao> 

Writing out half of Hamilton's equations for the full Hamiltonian of Eq. (23) 
we have 

dSt; 



dt 



i _1 rl , *m 

c 2 V 2 fi 



(31) 



relating the true velocity and pseudo velocity. The expressions for accel- 
erations in Sec. 4 become valid again if the v,- are interpreted as pseudo 
velocities. Operationally, this means that we have to transform the initial 
velocities to pseudo velocities by solving Eq. (31) for the v;, carry out the 
integration in terms of the pseudo velocities and then transform back to 
true velocities to output results. 

In our implementation, the relativistic corrections consume less than 
5% of the computing time. 



6. SAMPLE INTEGRATIONS 

We have implemented the individual time-step algorithm for nine planets. 
After some experimentation, we picked timestep ratios of 1 : 2 : 2 : 4 : 
8 : 8 : 64 : 64 : 256 for Mercury, Venus,. . ., Pluto, which makes the lon- 
gitudes roughly equally accurate for the planets. Accuracy of ~ 1 arcsec 
per century (or about 1 radian in 20 Myr) requires the smallest timestep 
to be about one week, assuming a warmup is used at the start. Symplec- 
tic interpolation was used as described in Sec. 4 to synchronize the Kepler 
clocks before advancing the interaction clocks. The computer time required 
is quite implementation and compiler dependent, but our code takes about 
15 sec for each kyr of integration on a 50 MHz Sparc-10; if all the timesteps 
are reduced to equal Mercury's our code takes about 35 sec per kyr. The 
latter case corresponds roughly to the Sussman & Wisdom (1992) inte- 
gration after refinements by WHT — the main difference is that the other 
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authors use a more approximate form for the post-Newtonian part. At the 
1 arcsec per century level the new algorithm is about an order of magni- 
tude faster than the 12th order symmetric multistep integrator in Quinlan 
& Tremaine (1990), but of course for high-enough accuracy requirements 
the latter would be more efficient (in the absence of roundoff errors). 

Here we check our implementation, with and without individual time- 
steps, against the 12th order symmetric multistep integrator. For the sym- 
plectic integrator we set Mercury's timestep to be 7-^ days; with the other 
timesteps in the ratios above a single cycle then spans 1800 days. For the 
warmup we first integrated back 5000 yr with all the timesteps reduced 
32-fold, while gradually reducing the interaction strengths to zero, then 
integrated forward to the starting point while gradually reviving the inter- 
actions again. For the multistep integrator we used a ^ day timestep, which 
should be much more accurate than the symplectic integrator. In Fig. 1 the 
symplectic integration has individual time-steps; in Fig. 2 the symplectic 
timesteps are all reduced to equal Mercury's. 

It is straightforward to estimate the time saved by integrating with 
individual timesteps. Suppose that advancing the Kepler clock K{ takes 
computer time AtK, and that advancing the interaction clock J; takes time 
(N—i)Ati (the number of interactions to be calculated equals the number of 
planets outside planet i). Then the ratio r of computer time to integration 
time is 



N 

E 



At K ,„ „Atj 



(N-i)^ 



T; 



(32) 



If a common timestep is used this expression simplifies to 



At K N(N-l)Atj , , 

r c = N— ~ + ^^ — . (33) 

r It 

For the timestep ratios we have used the reduction r/r c in computer time 
from using individual timesteps will lie between 0.28 (if Atx Ati) and 
0.46 (AtK "C At/); we actually obtained a reduction of 0.43. 

We offer no simple prescription for choosing the timesteps r;. The 
timesteps used in the sample integration were chosen by trial and error. 
Artificial stepsize resonances (Wisdom & Holman 1992) will be more nu- 
merous if there are individual timesteps, so comparison of two integrations 
with different timestep sets is always prudent. 

We have not addressed errors arising through roundoff, which can be 
significant in long integrations. The usual methods for controlling roundoff 
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10 4 10 5 10 6 

Time In yr 

1. Plotted here is the maximum error in mean anomaly M up to time T, against 
T, for leapfrog with individual timesteps. Mercury's timestep is 7^ days, and 
the other timesteps are larger in the ratios 1 : 2 : 2 : 4 : 8 : 8 : 64 : 64 : 256. 
The integration includes effects from general relativity and the Moon but does 
not integrate the lunar orbit. The faster-than-linear error growth for Mercury 
(noticeable for all the inner planets in Fig. 2) is presumably due to roundoff 
error. 

error, based in part on carrying out selected additions in quadruple precision 
(Quinn & Tremaine 1990), are not effective in MVS integrators, where most 
of the roundoff arises during the repeated conversion between Cartesian 
coordinates and Kepler elements. 
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Time in yr 

2. Like Fig. 1, but with all the planets having a timestep of 7^ days. Note 
that the vertical scale has been shrunk to accommodate the much larger range of 
accuracies. 
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7. VARIANTS 

Many variants are possible on the integrators we have discussed. A trivial 
and well-known example is to interchange the roles of the K and I clocks in 
the leapfrog cycle (10). For the two-timestep integrator (15), an alternative 
with similar properties is: 

(Advance K\ by ^r) 
( Advance I by r ) 

( Advance K\ by r and K2 by 2r ) (34) 
( Advance I by r ) 
(Advance K\ by ^r) 

In contrast to simple leapfrog, there is no time-reversible algorithm of equal 
simplicity to (15) or (34) that begins by advancing the clock I. 

The methods described here can be used to integrate other Hamilto- 
nians of the form (11) that satisfy the Poisson bracket relations (12). Our 
methods can even be applied to systems that do not satisfy the relations 
(12b): simply make the assignments i?; n t, 1 <— -ffint, #int,i <— 0, i = 2. .N . 

As an example of an application in another physical context, we con- 
sider the integrator for use in molecular dynamics that is described by Skeel 
k Biesiadecki (1994). They wish to follow the motion of a particle under 
the actions of forces F\ and F2, where F\ varies more rapidly than F2 and 
the timestep for evaluating F2 is M times the timestep for F\ (for simplic- 
ity we examine only the case M = 2). We can derive their algorithm from 
(15) by re-defining the clocks K{ and I: when K{ is running the particle's 
coordinates stay fixed while its momentum changes according to the force 
Fi, and while I is running its position changes at fixed momentum. For a 
particle of unit mass and phase-space coordinates (x,v) we then have 

( Increment v by \tF\ + ri<2 ) 
( Increment x by tv ) 

( Increment v by tF\ ) (35) 
( Increment x by tv ) 
( Increment v by \tF\ + ri<2 ) 
which is the integrator given by Skeel k Biesiadecki. 

We thank Pierre- Vincent Koseleff and Ken Wilson for discussions, and 
Mark Bartelt for advice on code optimization. This research was supported 
by NSERC. 
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